Here I study how biomass density changes across treatments in the PatchSizePilot. In particular, I’m studying how biomass density changes across:
culture_info = read.csv(here("data", "PatchSizePilot_culture_info.csv"), header = TRUE)
datatable(culture_info[,1:10],
rownames = FALSE,
options = list(scrollX = TRUE),
filter = list(position = 'top',
clear = FALSE))
Why did system number 40 disappear????
### --- IMPORT --- ###
load(here("data", "population", "t0.RData")); t0 = pop_output
load(here("data", "population", "t1.RData")); t1 = pop_output
load(here("data", "population", "t2.RData")); t2 = pop_output
load(here("data", "population", "t3.RData")); t3 = pop_output
load(here("data", "population", "t4.RData")); t4 = pop_output
load(here("data", "population", "t5.RData")); t5 = pop_output
load(here("data", "population", "t6.RData")); t6 = pop_output
load(here("data", "population", "t7.RData")); t7 = pop_output
rm(pop_output)
### --- TIDY --- ###
#Column: time
t0$time = NA
t1$time = NA
#Column: replicate_video
t0$replicate_video = 1:12 #In t1 I took 12 videos of a single
t1$replicate_video = 1 #In t1 I took only 1 video/culture
t2$replicate_video = 1 #In t2 I took only 1 video/culture
t3$replicate_video = 1 #In t3 I took only 1 video/culture
t4$replicate_video = 1 #In t4 I took only 1 video/culture
t5$replicate_video = 1 #In t5 I took only 1 video/culture
t6 = t6 %>%
rename(replicate_video = replicate)
t7 = t7 %>%
rename(replicate_video = replicate)
#Create an elongated version of t0 so that each of the 110 cultures can have 12 video replicates at t0.
elongating_t0 = NULL
for (video in 1:nrow(t0)){
for (ID in 1:nrow(culture_info)) {
elongating_t0 = rbind(elongating_t0, t0[video,])
}
}
ID_vector = rep(1:nrow(culture_info),
times = nrow(t0))
elongating_t0$culture_ID = ID_vector
#Merge previous data-sets
t0 = merge(culture_info,elongating_t0, by="culture_ID")
t1 = merge(culture_info,t1, by = "culture_ID")
t2 = merge(culture_info,t2, by = "culture_ID")
t3 = merge(culture_info,t3, by = "culture_ID")
t4 = merge(culture_info,t4, by = "culture_ID")
t5 = merge(culture_info,t5, by = "culture_ID")
t6 = merge(culture_info,t6, by = "culture_ID")
t7 = merge(culture_info,t7, by = "culture_ID")
ds_biomass = rbind(t0, t1, t2, t3, t4, t5, t6, t7)
rm(elongating_t0, t0, t1, t2, t3, t4, t5, t6, t7)
ds_biomass$time_point[ds_biomass$time_point=="t0"] = 0
ds_biomass$time_point[ds_biomass$time_point=="t1"] = 1
ds_biomass$time_point[ds_biomass$time_point=="t2"] = 2
ds_biomass$time_point[ds_biomass$time_point=="t3"] = 3
ds_biomass$time_point[ds_biomass$time_point=="t4"] = 4
ds_biomass$time_point[ds_biomass$time_point=="t5"] = 5
ds_biomass$time_point[ds_biomass$time_point=="t6"] = 6
ds_biomass$time_point[ds_biomass$time_point=="t7"] = 7
ds_biomass$time_point = as.character(ds_biomass$time_point)
ds_biomass$day = NA
ds_biomass$day[ds_biomass$time_point== 0] = 0
ds_biomass$day[ds_biomass$time_point== 1] = 4
ds_biomass$day[ds_biomass$time_point== 2] = 8
ds_biomass$day[ds_biomass$time_point== 3] = 12
ds_biomass$day[ds_biomass$time_point== 4] = 16
ds_biomass$day[ds_biomass$time_point== 5] = 20
ds_biomass$day[ds_biomass$time_point== 6] = 24
ds_biomass$day[ds_biomass$time_point== 7] = 28
#Column: eco_metaeco_type
ds_biomass$eco_metaeco_type = factor(ds_biomass$eco_metaeco_type,
levels = c('S',
'S (S_S)',
'S (S_L)',
'M',
'M (M_M)',
'L',
'L (L_L)',
'L (S_L)'))
ecosystems_to_take_off = 60 #Culture number 60 because it was spilled
ds_biomass = ds_biomass %>%
filter(! culture_ID %in% ecosystems_to_take_off)
ds_for_evaporation = ds_biomass
ds_biomass = ds_biomass %>%
select(culture_ID,
patch_size,
disturbance,
metaecosystem_type,
bioarea_per_volume,
replicate_video,
time_point,
day,
metaecosystem,
system_nr,
eco_metaeco_type) %>%
relocate(culture_ID,
system_nr,
disturbance,
time_point,
day,
patch_size,
metaecosystem,
metaecosystem_type,
eco_metaeco_type,
replicate_video,
bioarea_per_volume)
datatable(ds_biomass,
rownames = FALSE,
options = list(scrollX = TRUE),
filter = list(position = 'top',
clear = FALSE))
ds_regional = ds_biomass %>%
filter(metaecosystem == "yes") %>%
group_by(culture_ID,
system_nr,
disturbance,
day,
time_point,
patch_size,
metaecosystem_type) %>%
summarise(patch_mean_bioarea_across_videos = mean(bioarea_per_volume)) %>%
group_by(system_nr, disturbance, day, time_point, metaecosystem_type) %>%
summarise(regional_mean_bioarea = mean(patch_mean_bioarea_across_videos))
metaecosystems_to_take_off = 40 #System 40 was the system of culture 60 that I spilled
ds_regional = ds_regional %>%
filter(! system_nr %in% metaecosystems_to_take_off)
datatable(ds_regional,
rownames = FALSE,
options = list(scrollX = TRUE),
filter = list(position = 'top',
clear = FALSE))
We want to understand if the regional biomass produced by an ecosystem with a small and a large patch (metaecosystem_type = S_L) is lower than the regional biomass produced by an ecosystem with two medium patches (metaecosystem_type = M_M).
Let’s start by plotting the single ecosystems to see that everything is fine. To make the patterns clear let’s plot the low disturbance and high disturbance in two different figures. We first plot the single meta-ecosystems and then their box plots.
ds_regional %>%
filter ( disturbance == "low") %>%
filter (metaecosystem_type == "S_L" |
metaecosystem_type == "M_M") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = metaecosystem_type)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = low",
fill = "System nr",
color = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_linetype_discrete(labels = c("medium-medium",
"small-large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional %>%
filter ( disturbance == "high") %>%
filter (metaecosystem_type == "S_L" | metaecosystem_type == "M_M") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = metaecosystem_type)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = high",
fill = "System nr",
color = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_linetype_discrete(labels = c("medium-medium",
"small-large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
p_regional_low_mean = ds_regional %>%
filter(disturbance == "low") %>%
filter (metaecosystem_type == "S_L" | metaecosystem_type == "M_M") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = interaction(day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = low",
color='',
fill='') +
scale_fill_discrete(labels = c("medium-medium",
"small-large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
p_regional_low_mean
ds_regional %>%
filter(disturbance == "high") %>%
filter (metaecosystem_type == "S_L" | metaecosystem_type == "M_M") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = interaction (day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = high",
color='',
fill='') +
scale_fill_discrete(labels = c("medium-medium", "small-large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
We can see that the regional biomass was higher for meta-ecosystems with the same size, regardless of disturbance. As both disturbance levels showed the same pattern, I would just keep the one with low disturbance in the publication, as it shows a stronger pattern.
Tidy
First of all, let’s modify the data set including the regional biomass of our meta-ecosystems. In this data set, we want to have the regional biomass of the meta-ecosystems (averaged first across videos and then across patches) in which we:
Include only the meta-ecosystems in which patches had both medium size (metaecosystem_type = M_M) and meta-ecosystems in which patches had one a small size and the other large size (metaecosystem_type = S_L).
Take off the first two point (day 0 and day = 4). This is because the first perturbation happened only at day 5.
ds_regional_MM_SL_t2t7 = ds_regional %>%
filter (metaecosystem_type == "M_M" | metaecosystem_type == "S_L",
time_point >= 2)
Model selection
Let’s start from the largest mixed effect model.
full_model = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
metaecosystem_type : disturbance +
(metaecosystem_type | day) +
(disturbance | day) +
(metaecosystem_type : disturbance | day) +
(1 | system_nr) ,
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
Should we keep M * D?
no_MD = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
(metaecosystem_type | day) +
(disturbance | day) +
(1 | system_nr) ,
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
anova(full_model, no_MD)
## Data: ds_regional_MM_SL_t2t7
## Models:
## no_MD: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
## full_model: regional_mean_bioarea ~ metaecosystem_type + disturbance + metaecosystem_type:disturbance + (metaecosystem_type | day) + (disturbance | day) + (metaecosystem_type:disturbance | day) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_MD 11 1861.7 1892.4 -919.86 1839.7
## full_model 27 1891.1 1966.4 -918.56 1837.1 2.5891 16 0.9999
No.
Should we keep the random slope of (M | day)?
no_M_day_slope = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
(disturbance | day) +
(1 | system_nr) ,
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
anova(no_MD, no_M_day_slope)
## Data: ds_regional_MM_SL_t2t7
## Models:
## no_M_day_slope: regional_mean_bioarea ~ metaecosystem_type + disturbance + (disturbance | day) + (1 | system_nr)
## no_MD: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_M_day_slope 8 1871.0 1893.2 -927.48 1855.0
## no_MD 11 1861.7 1892.4 -919.86 1839.7 15.239 3 0.001623 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yes.
Should we keep the correlation between intercept and slope of (M | day)?
no_M_day_correlation = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
(metaecosystem_type || day) +
(disturbance | day) +
(1 | system_nr) ,
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
anova(no_MD, no_M_day_correlation)
## Data: ds_regional_MM_SL_t2t7
## Models:
## no_MD: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
## no_M_day_correlation: regional_mean_bioarea ~ metaecosystem_type + disturbance + ((1 | day) + (0 + metaecosystem_type | day)) + (disturbance | day) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_MD 11 1861.7 1892.4 -919.86 1839.7
## no_M_day_correlation 12 1863.7 1897.2 -919.86 1839.7 0 1 0.9996
Yes.
Should we keep the random slope of (D| day)?
no_D_day_slope = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
(metaecosystem_type | day) +
(1 | system_nr) ,
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
anova(no_MD, no_M_day_slope)
## Data: ds_regional_MM_SL_t2t7
## Models:
## no_M_day_slope: regional_mean_bioarea ~ metaecosystem_type + disturbance + (disturbance | day) + (1 | system_nr)
## no_MD: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_M_day_slope 8 1871.0 1893.2 -927.48 1855.0
## no_MD 11 1861.7 1892.4 -919.86 1839.7 15.239 3 0.001623 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yes.
Should we keep the correlation between intercept and slope of (D | day)?
no_D_day_correlation = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
(metaecosystem_type | day) +
(disturbance || day) +
(1 | system_nr) ,
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
anova(no_MD, no_D_day_correlation)
## Data: ds_regional_MM_SL_t2t7
## Models:
## no_MD: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + (disturbance | day) + (1 | system_nr)
## no_D_day_correlation: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + ((1 | day) + (0 + disturbance | day)) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_MD 11 1861.7 1892.4 -919.86 1839.7
## no_D_day_correlation 12 1863.7 1897.2 -919.86 1839.7 0 1 0.9996
No.
Should we keep (1 | system_nr)?
no_system_nr = lmer(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
(metaecosystem_type | day) +
(disturbance || day),
data = ds_regional_MM_SL_t2t7,
REML = FALSE)
anova(no_D_day_correlation, no_system_nr)
## Data: ds_regional_MM_SL_t2t7
## Models:
## no_system_nr: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + ((1 | day) + (0 + disturbance | day))
## no_D_day_correlation: regional_mean_bioarea ~ metaecosystem_type + disturbance + (metaecosystem_type | day) + ((1 | day) + (0 + disturbance | day)) + (1 | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_system_nr 11 1869.2 1899.9 -923.60 1847.2
## no_D_day_correlation 12 1863.7 1897.2 -919.86 1839.7 7.4907 1 0.006202
##
## no_system_nr
## no_D_day_correlation **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yes.
Best model
Therefore, the best model is
\[ Regional \: bioarea = M + D + (M | t) + (D || t) \]
This model is the best model when looking at all time points coming after the first disturbance event (t2->t7). Assuming that this model holds for also other sections of the time series, the r squared of the model and of meta-ecosystem type is as follows.
#Create a table with all the models in which time is a random effect.
### --- INITIALISE TABLE --- ###
columns = c("model", "time_point", "AIC", "BIC", "R2_mixed", "R2_fixed", "R2_mixed_M", "R2_fixed_M")
random_time_table = data.frame(matrix(ncol = length(columns),
nrow = 0))
colnames(random_time_table) = columns
### --- M + D + (M | t) + (D || t) --- ###
for (last_point in 3:7) {
full_model = lmer(regional_mean_bioarea ~
disturbance +
metaecosystem_type +
(metaecosystem_type | day) +
(1 | system_nr),
data = filter(ds_regional_MM_SL_t2t7,
time_point <= last_point),
REML = FALSE,
control = lmerControl(optimizer ="Nelder_Mead"))
null_model = lm(regional_mean_bioarea ~
1 ,
data = filter(ds_regional_MM_SL_t2t7, time_point <= last_point))
metaeco_null = lmer(regional_mean_bioarea ~
disturbance +
(1 | day) +
(1 | system_nr),
data = filter(ds_regional_MM_SL_t2t7,
time_point <= last_point),
REML = FALSE,
control = lmerControl(optimizer ="Nelder_Mead"))
random_time_table = update_all_models_table("M+D+(M|t)+(D||t)",
random_time_table,
full_model,
null_model,
metaeco_null,
"mixed")
}
datatable(random_time_table,
rownames = FALSE,
options = list(pageLength = 100,
scrollX = TRUE,
autoWidth = TRUE,
columnDefs = list(list(targets=c(0),visible=TRUE, width='160'),
list(targets=c(1), visible=TRUE, width='10'),
list(targets=c(2), visible=TRUE, width='10'),
list(targets=c(3), visible=TRUE, width='10'),
list(targets=c(4), visible=TRUE, width='10'),
list(targets=c(5), visible=TRUE, width='10'),
list(targets=c(6), visible=TRUE, width='10'),
list(targets=c(7), visible=TRUE, width='10'),
list(targets='_all', visible=FALSE))),
caption = "
M = Meta-ecosystem type,
D = disturbance,
t = time,
(M | t) = random effect of time on the intercept and slope of M,
(D || t) = random effect of time on the intercept and slope of D,
|| = no correlation between intercept and slope,
| = correlation between intercept and slope,
R2_mixed = r squared of the model,
R2_fixed = r squared of the model when considering only fixed effects,
R2_mixed_M = r squared of meta-ecosystem type,
R2_fixed_M = r squared of meta-ecosystem type when considerin only fixed effects")
How is it possible that the mixed effect of M is higher than its fixed effects?
Here we want to see whether we can include time as a fixed effect by transforming the bioarea into its logarithm with base 10. As we are not considering the first two time points because they were before the first disturbance (we want to see the effects of meta-ecosystem type under disturbance), we are lucky that there is a better chance that the biomass will look like goes down in a linear fashion.
Linearity of regional bioarea ~ time
ds_regional %>%
filter(time_point >= 2) %>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = day)) +
geom_boxplot() +
labs(title = "Without log transformation",
x = "Day",
y = "Regional bioarea (something/µl)")
Let’s check how linear the relationship is.
linear_model = lm(regional_mean_bioarea ~
day,
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"))
par(mfrow=c(2,3))
plot(linear_model, which = 1:5)
I’m not 100% convinced, as the residuals vs fitted has a bit of a banana shape line.
Linearity of Log10(Regional bioarea +1) ~ time
ds_regional %>%
filter(time_point >= 2) %>%
ggplot(aes(x = day,
y = log(regional_mean_bioarea + 1),
group = day)) +
geom_boxplot() +
labs(title = "With log transformation",
x = "Day",
y = "Log (regional bioarea + 1) (something/µl)")
Let’s now check how linear these two relationships are.
log_linear_model = lm(log10(regional_mean_bioarea + 1) ~
day,
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"))
par(mfrow=c(2,3))
plot(log_linear_model, which = 1:5)
par(mfrow=c(1,1))
Way better. Especially the residuals vs fitted values plot doesn’t have a banana shape anymore. This seems to be a good model. Let’s then keep the log transformed bioarea.
Model selection
Let’ start from the full model.
full = lmer(log10(regional_mean_bioarea + 1) ~
day * metaecosystem_type * disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
Should we keep the correlation in (day | system_nr)?
no_correlation = lmer(log10(regional_mean_bioarea + 1) ~
day * metaecosystem_type * disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
anova(full, no_correlation)
## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## full: log10(regional_mean_bioarea + 1) ~ day * metaecosystem_type * disturbance + (day | system_nr)
## no_correlation: log10(regional_mean_bioarea + 1) ~ day * metaecosystem_type * disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## full 12 -168.78 -135.33 96.389 -192.78
## no_correlation 12 -168.78 -135.33 96.389 -192.78 0 0
Yes.
Should we keep t * M * D?
no_threeway = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
day : metaecosystem_type +
day : disturbance +
metaecosystem_type : disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = 'optimx',
optCtrl = list(method = 'L-BFGS-B')))
anova(full, no_threeway)
## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## no_threeway: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:metaecosystem_type + day:disturbance + metaecosystem_type:disturbance + (day | system_nr)
## full: log10(regional_mean_bioarea + 1) ~ day * metaecosystem_type * disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_threeway 11 -170.66 -140.00 96.330 -192.66
## full 12 -168.78 -135.33 96.389 -192.78 0.1182 1 0.7309
No.
Should we keep t * M?
no_TM = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
day : disturbance +
metaecosystem_type : disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
anova(no_threeway,no_TM)
## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## no_TM: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + metaecosystem_type:disturbance + (day | system_nr)
## no_threeway: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:metaecosystem_type + day:disturbance + metaecosystem_type:disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_TM 10 -172.62 -144.74 96.308 -192.62
## no_threeway 11 -170.66 -140.00 96.330 -192.66 0.0431 1 0.8356
No.
Should we keep t * D?
no_TD = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
metaecosystem_type : disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
anova(no_TM, no_TD)
## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## no_TD: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + metaecosystem_type:disturbance + (day | system_nr)
## no_TM: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + metaecosystem_type:disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_TD 9 -163.54 -138.46 90.771 -181.54
## no_TM 10 -172.62 -144.74 96.308 -192.62 11.074 1 0.0008756 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No.
Should we keep M * D?
no_MD = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
day : disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
anova(no_TM, no_MD)
## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## no_MD: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + (day | system_nr)
## no_TM: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + metaecosystem_type:disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_MD 9 -173.87 -148.78 95.932 -191.87
## no_TM 10 -172.62 -144.74 96.308 -192.62 0.7513 1 0.3861
No.
Should we keep the random effect of system nr on the time slopes (day | system_nr)?
no_random_slopes = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
day : disturbance +
(1 | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
anova(no_MD, no_random_slopes)
## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## no_random_slopes: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + (1 | system_nr)
## no_MD: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_random_slopes 7 -173.23 -153.71 93.613 -187.23
## no_MD 9 -173.87 -148.78 95.932 -191.87 4.6394 2 0.0983 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yes.
Best model
Therefore, our best model is:
\[ log_{10}(regional \: bioarea + 1) = t + M + D + t*D + (t| system \: nr) \]
This model is the best model when looking at all time points coming after the first disturbance event (t2->t7). Assuming that this model holds for also other sections of the time series, the r squared of the model and of meta-ecosystem type is as follows.
#Create a table in which the regional biomass has been log transformed.
### --- INITIALISE TABLE --- ###
columns = c("model", "time_point", "AIC", "BIC", "R2_mixed", "R2_fixed", "R2_mixed_M", "R2_fixed_M")
log_time_table = data.frame(matrix(ncol = length(columns), nrow = 0))
colnames(log_time_table) = columns
### --- POPULATE THE TABLE --- ###
for (last_point in 4:7) {
full_model = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
day : disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(time_point <= last_point) %>%
filter(metaecosystem_type == "M_M" |
metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
null_model = lm(regional_mean_bioarea ~
1 ,
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(time_point <= last_point) %>%
filter(metaecosystem_type == "M_M" |
metaecosystem_type == "S_L"))
metaeco_null_model = lmer(log10(regional_mean_bioarea + 1) ~
day +
disturbance +
day : disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(time_point <= last_point) %>%
filter(metaecosystem_type == "M_M" |
metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
log_time_table = update_all_models_table("t + M + D + t * M * D + (t || system_nr)",
log_time_table,
full_model,
null_model,
metaeco_null_model,
"mixed")
}
datatable(log_time_table,
rownames = FALSE,
options = list(pageLength = 100,
scrollX = TRUE,
autoWidth = TRUE,
columnDefs = list(list(targets=c(0),visible=TRUE, width='160'),
list(targets=c(1), visible=TRUE, width='10'),
list(targets=c(2), visible=TRUE, width='10'),
list(targets=c(3), visible=TRUE, width='10'),
list(targets=c(4), visible=TRUE, width='10'),
list(targets=c(5), visible=TRUE, width='10'),
list(targets=c(6), visible=TRUE, width='10'),
list(targets=c(7), visible=TRUE, width='10'),
list(targets='_all', visible=FALSE))),
caption = "
M = Meta-ecosystem type,
D = disturbance,
(1 | t) = random effect of time on the intercept,
(1 | ID) = random effect of meta-ecosystem ID on the intercept,
|| = no correlation between intercept and slope,
| = correlation between intercept and slope,
R2 = r squared of the whole model,
R2_fixed = fixed part of the mixed model,
mixed_R2 = r squared when considering both fixed and random effects (conditional r squared),
fixed_R2 = r squared when considering only the fixed effects (marginal r squared)")
Notice that I did not include the t3-t4 time series, as when running the model it gives me the following error:
Here we want to fit how biomass changes across time to a function. The biomass of our meta-ecosystems looks like this.
ds_regional %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L") %>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = day)) +
geom_boxplot() +
labs(x = "day", y = "Regional bioarea (something/microlitres)") +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
To fit these data, we need to produce (and parameterise) a function that resemble how the biomass first increases and then decreases.
Hank produced the following function:
\[biomass = a_4 * (day - a_5) * e^{a_1(day - a_5)}\]
If we parameterise the function and then plot, it looks as follows.
a1 = -0.1
a4 = 1200
a5 = -1
day = seq(0, 30, 0.01)
biomass = a4*(day-a5) * exp(a1*(day-a5))
plot(biomass ~ day)
Now, let’s find the best parameters (a1, a4, a5) that fit our data.
ds_regional_shrunk_type = ds_regional %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L")
model = nls(regional_mean_bioarea ~ a4 * (day-a5) * exp(a1 * (day-a5)),
start = list(a1 = -0.1, a4 = 1200, a5 = -1),
trace = T,
data = ds_regional_shrunk_type)
a1 = as.numeric(model$m$getPars()[1])
a4 = as.numeric(model$m$getPars()[2])
a5 = as.numeric(model$m$getPars()[3])
model$m$getPars()
## a1 a4 a5
## -0.1336182 1216.4268967 -1.6147345
And now let’s plot the function to see how it fits our data.
day = seq(0,30,1)
predicted = a4*(day-a5)*exp(a1*(day-a5))
data_fitted=data.frame(day=day,regional_mean_bioarea=predicted)
ds_regional_shrunk_type%>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = day)) +
geom_boxplot() +
labs(x = "Day", y = "Regional bioarea") +
geom_line(data=data_fitted,aes(x = day, y=regional_mean_bioarea),color="red", group = 1) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
Let’s then include our predictions into the data set.
ds_regional_predicted_shrunk_type = ds_regional %>%
mutate(predicted_from_time = a4*(day-a5)*exp(a1*(day-a5))) %>%
filter(metaecosystem_type == "S_L" | metaecosystem_type == "M_M")
Let’s now work using the fitted data that we found in the previous section.
Tidy
ds_regional_predicted_shrunk_type_n_day = ds_regional_predicted_shrunk_type %>%
filter(time_point >= 2)
Model selection
Let’s start from the largest mixed effect model.
full = lmer(log10(regional_mean_bioarea + 1) ~
predicted_from_time * metaecosystem_type * disturbance +
(predicted_from_time | system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
## Warning: Some predictor variables are on very different scales: consider
## rescaling
Should we keep the correlation in (day | system_nr)?
no_correlation = lmer(log10(regional_mean_bioarea + 1) ~
predicted_from_time * metaecosystem_type * disturbance +
(predicted_from_time || system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
## Warning: Some predictor variables are on very different scales: consider
## rescaling
anova(full, no_correlation)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## no_correlation: log10(regional_mean_bioarea + 1) ~ predicted_from_time * metaecosystem_type * disturbance + ((1 | system_nr) + (0 + predicted_from_time | system_nr))
## full: log10(regional_mean_bioarea + 1) ~ predicted_from_time * metaecosystem_type * disturbance + (predicted_from_time | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_correlation 11 -158.51 -127.85 90.257 -180.51
## full 12 -161.16 -127.71 92.581 -185.16 4.6486 1 0.03108 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Keep full
Yes.
Should we keep t * M * D?
no_TMD = lmer(log10(regional_mean_bioarea + 1) ~
predicted_from_time +
metaecosystem_type +
disturbance +
predicted_from_time : metaecosystem_type +
predicted_from_time : disturbance +
metaecosystem_type : disturbance +
(predicted_from_time | system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
## Warning: Some predictor variables are on very different scales: consider
## rescaling
anova(full, no_TMD)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## no_TMD: log10(regional_mean_bioarea + 1) ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time:metaecosystem_type + predicted_from_time:disturbance + metaecosystem_type:disturbance + (predicted_from_time | system_nr)
## full: log10(regional_mean_bioarea + 1) ~ predicted_from_time * metaecosystem_type * disturbance + (predicted_from_time | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_TMD 11 -163.07 -132.41 92.538 -185.07
## full 12 -161.16 -127.71 92.581 -185.16 0.0874 1 0.7675
#Keep no_TMD
No.
Should we keep t * M?
no_TM = lmer(log10(regional_mean_bioarea + 1) ~
predicted_from_time +
metaecosystem_type +
disturbance +
predicted_from_time : disturbance +
metaecosystem_type : disturbance +
(predicted_from_time | system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
## Warning: Some predictor variables are on very different scales: consider
## rescaling
anova(no_TMD, no_TM)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## no_TM: log10(regional_mean_bioarea + 1) ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time:disturbance + metaecosystem_type:disturbance + (predicted_from_time | system_nr)
## no_TMD: log10(regional_mean_bioarea + 1) ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time:metaecosystem_type + predicted_from_time:disturbance + metaecosystem_type:disturbance + (predicted_from_time | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_TM 10 -165.06 -137.18 92.528 -185.06
## no_TMD 11 -163.07 -132.41 92.538 -185.07 0.0194 1 0.8893
#Keep no_TM
No.
Should we keep t * D?
no_TD = lmer(log10(regional_mean_bioarea + 1) ~
predicted_from_time +
metaecosystem_type +
disturbance +
metaecosystem_type : disturbance +
(predicted_from_time | system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
anova(no_TM, no_TD)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## no_TD: log10(regional_mean_bioarea + 1) ~ predicted_from_time + metaecosystem_type + disturbance + metaecosystem_type:disturbance + (predicted_from_time | system_nr)
## no_TM: log10(regional_mean_bioarea + 1) ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time:disturbance + metaecosystem_type:disturbance + (predicted_from_time | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_TD 9 -156.84 -131.75 87.420 -174.84
## no_TM 10 -165.06 -137.18 92.528 -185.06 10.216 1 0.001392 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Keep no_TM
Yes.
Should we keep M * D?
no_MD = lmer(log10(regional_mean_bioarea + 1) ~
predicted_from_time +
metaecosystem_type +
disturbance +
predicted_from_time : disturbance +
(predicted_from_time | system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
## Warning: Some predictor variables are on very different scales: consider
## rescaling
anova(no_TM, no_MD)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## no_MD: log10(regional_mean_bioarea + 1) ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time:disturbance + (predicted_from_time | system_nr)
## no_TM: log10(regional_mean_bioarea + 1) ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time:disturbance + metaecosystem_type:disturbance + (predicted_from_time | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_MD 9 -166.26 -141.17 92.129 -184.26
## no_TM 10 -165.06 -137.18 92.528 -185.06 0.7978 1 0.3717
#Keep no_MD
No.
Should we keep the random effect of system nr on the time slopes (day | system_nr)?
no_random_slopes = lmer(log10(regional_mean_bioarea + 1) ~
predicted_from_time +
metaecosystem_type +
disturbance +
predicted_from_time : disturbance +
(1 | system_nr),
data = ds_regional_predicted_shrunk_type_n_day,
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
## Warning: Some predictor variables are on very different scales: consider
## rescaling
anova(no_MD, no_random_slopes)
## Data: ds_regional_predicted_shrunk_type_n_day
## Models:
## no_random_slopes: log10(regional_mean_bioarea + 1) ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time:disturbance + (1 | system_nr)
## no_MD: log10(regional_mean_bioarea + 1) ~ predicted_from_time + metaecosystem_type + disturbance + predicted_from_time:disturbance + (predicted_from_time | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_random_slopes 7 -165.49 -145.98 89.746 -179.49
## no_MD 9 -166.26 -141.17 92.129 -184.26 4.7662 2 0.09226 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Keep no_MD
Yes.
Best model
Therefore, our best model is:
\[ log_{10}(regional \: bioarea + 1) = t + M + D + t*D + (t| system \: nr) \]
This model is the best model when looking at all time points coming after the first disturbance event (t2->t7). Assuming that this model holds for also other sections of the time series, the r squared of the model and of meta-ecosystem type is as follows.
#Create a table in which time is a fixed effect.
### --- INITIALISE TABLE --- ###
columns = c("model", "time_point", "AIC", "BIC", "R2_mixed", "R2_fixed", "R2_mixed_M", "R2_fixed_M")
fitted_time_table = data.frame(matrix(ncol = length(columns), nrow = 0))
colnames(fitted_time_table) = columns
### --- POPULATE TABLE --- ###
for (last_point in 4:7) {
full_model = lmer(log10(regional_mean_bioarea + 1) ~
predicted_from_time +
metaecosystem_type +
disturbance +
predicted_from_time : disturbance +
(predicted_from_time | system_nr),
data = ds_regional_predicted_shrunk_type_n_day %>%
filter(time_point <= last_point),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
null_model = lm(regional_mean_bioarea ~
1 ,
data = ds_regional_predicted_shrunk_type_n_day %>%
filter(time_point <= last_point))
metaeco_null_model = lmer(log10(regional_mean_bioarea + 1) ~
predicted_from_time +
disturbance +
predicted_from_time : disturbance +
(predicted_from_time | system_nr),
data = ds_regional_predicted_shrunk_type_n_day %>%
filter(time_point <= last_point),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
fitted_time_table = update_all_models_table("Tp + M + D + Tp * D + (Tp | system_nr)",
fitted_time_table,
full_model,
null_model,
metaeco_null_model,
"mixed")
}
datatable(fitted_time_table,
rownames = FALSE,
options = list(pageLength = 100,
scrollX = TRUE,
autoWidth = TRUE,
columnDefs = list(list(targets=c(0),visible=TRUE, width='160'),
list(targets=c(1), visible=TRUE, width='10'),
list(targets=c(2), visible=TRUE, width='10'),
list(targets=c(3), visible=TRUE, width='10'),
list(targets=c(4), visible=TRUE, width='10'),
list(targets=c(5), visible=TRUE, width='10'),
list(targets=c(6), visible=TRUE, width='10'),
list(targets=c(7), visible=TRUE, width='10'),
list(targets='_all', visible=FALSE))),
caption = "
Tp = predicted from time,
M = Meta-ecosystem type,
D = disturbance,
(1 | t) = random effect of time on the intercept,
(1 | ID) = random effect of meta-ecosystem ID on the intercept,
|| = no correlation between intercept and slope,
| = correlation between intercept and slope,
R2 = r squared of the whole model,
R2_fixed = fixed part of the mixed model,
mixed_R2 = r squared when considering both fixed and random effects (conditional r squared),
fixed_R2 = r squared when considering only the fixed effects (marginal r squared)")
#Table with all the different models
full_table = rbind(random_time_table, fitted_time_table, log_time_table)
datatable(full_table,
rownames = FALSE,
options = list(pageLength = 100,
scrollX = TRUE,
autoWidth = TRUE,
columnDefs = list(list(targets=c(0),visible=TRUE, width='160'),
list(targets=c(1), visible=TRUE, width='10'),
list(targets=c(2), visible=TRUE, width='10'),
list(targets=c(3), visible=TRUE, width='10'),
list(targets=c(4), visible=TRUE, width='10'),
list(targets=c(5), visible=TRUE, width='10'),
list(targets=c(6), visible=TRUE, width='10'),
list(targets=c(7), visible=TRUE, width='10'),
list(targets='_all', visible=FALSE))),
caption = "
M = Meta-ecosystem type,
D = disturbance,
(1 | t) = random effect of time on the intercept,
(1 | ID) = random effect of meta-ecosystem ID on the intercept,
|| = no correlation between intercept and slope,
| = correlation between intercept and slope,
R2 = r squared of the whole model,
R2_fixed = fixed part of the mixed model,
mixed_R2 = r squared when considering both fixed and random effects (conditional r squared),
fixed_R2 = r squared when considering only the fixed effects (marginal r squared)")
Model selection
I’ll do model selection only on time point number 3 (however, I have done it also with time point 4,5,6,7 and they all give me the same result). Let’s start from the full model.
I can’t construct it from a mixed model because the following error pops up:
full = lm(regional_mean_bioarea ~
metaecosystem_type +
disturbance +
disturbance : metaecosystem_type,
data = ds_regional %>%
filter(time_point == 3) %>%
filter(metaecosystem_type == "M_M" |
metaecosystem_type == "S_L"))
Should we keep D * M?
no_MD = lm(regional_mean_bioarea ~
metaecosystem_type +
disturbance,
data = ds_regional %>%
filter(time_point == 3) %>%
filter(metaecosystem_type == "M_M" |
metaecosystem_type == "S_L"))
anova(full, no_MD)
## Analysis of Variance Table
##
## Model 1: regional_mean_bioarea ~ metaecosystem_type + disturbance + disturbance:metaecosystem_type
## Model 2: regional_mean_bioarea ~ metaecosystem_type + disturbance
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 16 3433038
## 2 17 3732239 -1 -299201 1.3945 0.2549
AIC(full, no_MD)
## df AIC
## full 5 307.8220
## no_MD 4 307.4933
No.
Best model
Therefore, our best model is
\[ Regional \: Bioarea = M + D \]
columns = c("time_point", "R2", "R2_M")
single_points = matrix(ncol = length(columns),
nrow = 7)
single_points = as.data.frame(single_points)
colnames(single_points) = columns
for (t in 1:7) {
full_model = lm(regional_mean_bioarea ~
disturbance +
metaecosystem_type,
data = ds_regional %>%
filter(metaecosystem_type == "M_M" |
metaecosystem_type == "S_L") %>%
filter(time_point == t))
no_M_model = lm(regional_mean_bioarea ~
disturbance,
data = ds_regional %>%
filter(metaecosystem_type == "M_M" |
metaecosystem_type == "S_L") %>%
filter(time_point == t))
R2_full_model = summary(full_model)$adj.r.squared
R2_no_M_model = summary(no_M_model)$adj.r.squared
R2_M = R2_full_model - R2_no_M_model
single_points$time_point[t] = t
single_points$R2[t] = R2_full_model
single_points$R2_M[t] = R2_M
}
single_points = round(single_points, digits = 2)
single_points = single_points[2:nrow(single_points),]
datatable(single_points,
rownames = FALSE,
colnames = c("Time point", "R2 model", "R2 meta-ecosystem type"))
Here I test whether the bioarea of small-large we see is due to meta-ecosystem dynamics or simply to their area.
To make the data-set, I’m going through the following steps:
isolated_S = ds_biomass %>%
filter(eco_metaeco_type == "S") %>%
group_by(system_nr, disturbance, time_point, day, patch_size) %>%
summarise(bioarea_per_volume_across_videos = mean(bioarea_per_volume))
isolated_L = ds_biomass %>%
filter(eco_metaeco_type == "L") %>%
group_by(system_nr, disturbance, time_point, day, patch_size) %>%
summarise(bioarea_per_volume_across_videos = mean(bioarea_per_volume))
### Check that they have the same number of patches, then check which one is the oen missing
length(unique(isolated_S$system_nr))
length(unique(isolated_L$system_nr))
unique(isolated_S$system_nr)
unique(isolated_L$system_nr)
#Take off one of the high disturbance small patches (system nr 50) because isolated small misses one at high disturbance
isolated_L = isolated_L%>%
filter(!system_nr == 50)
n_isolated_patches = length(unique(isolated_L$system_nr))
n_time_points = 8
number_for_pairing = rep( c( 1:n_isolated_patches), each = n_time_points)
number_for_pairing = as.data.frame(number_for_pairing)
colnames(number_for_pairing) = "number_for_pairing"
isolated_S = cbind(isolated_S, number_for_pairing)
isolated_L = cbind(isolated_L, number_for_pairing)
isolated_S_and_L = rbind(isolated_S, isolated_L)
SL_from_isolated = isolated_S_and_L %>%
group_by(disturbance, time_point, day, number_for_pairing) %>%
summarise(regional_mean_bioarea = mean(bioarea_per_volume_across_videos)) %>%
mutate(metaecosystem_type = "S_L_from_isolated") %>%
mutate(system_nr = number_for_pairing)
ds_regional_with_SL_from_isolated = rbind(SL_from_isolated, ds_regional)
ds_regional_with_SL_from_isolated %>%
filter(disturbance == "low") %>%
filter(metaecosystem_type == "S_L" | metaecosystem_type == "S_L_from_isolated") %>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = interaction(day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(title = "Disturbance = low",
x = "Day",
y = "Regional bioarea (something/microlitre)",
fill = "") +
scale_fill_discrete(labels = c("small-large", "isolated small & \n isolated large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional_with_SL_from_isolated %>%
filter(disturbance == "high") %>%
filter(metaecosystem_type == "S_L" | metaecosystem_type == "S_L_from_isolated") %>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = interaction(day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(title = "Disturbance = high",
x = "Day",
y = "Regional bioarea (something/microlitre)",
fill = "") +
scale_fill_discrete(labels = c("small-large", "isolated small & \n isolated large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
mixed_model = lmer(log10(regional_mean_bioarea +1) ~
day * metaecosystem_type * disturbance +
(day | system_nr),
data = ds_regional_with_SL_from_isolated %>%
filter(metaecosystem_type == "S_L" |
metaecosystem_type == "S_L_from_isolated") %>%
filter(time_point >= 2),
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
null_model = lmer(log10(regional_mean_bioarea +1) ~
day * disturbance +
(day | system_nr),
data = ds_regional_with_SL_from_isolated %>%
filter(metaecosystem_type == "S_L" |
metaecosystem_type == "S_L_from_isolated") %>%
filter(time_point >= 2),
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
anova(mixed_model, null_model)
## Data: ds_regional_with_SL_from_isolated %>% filter(metaecosystem_type == ...
## Models:
## null_model: log10(regional_mean_bioarea + 1) ~ day * disturbance + (day | system_nr)
## mixed_model: log10(regional_mean_bioarea + 1) ~ day * metaecosystem_type * disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## null_model 8 -154.45 -132.56 85.224 -170.45
## mixed_model 12 -147.96 -115.12 85.980 -171.96 1.5111 4 0.8247
ds_regional %>%
filter(!metaecosystem_type == "S_L") %>%
filter ( disturbance == "low") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = metaecosystem_type)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = low",
fill = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_colour_continuous(guide = "none") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large-large",
"medium-medium",
"small-small")) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional %>%
filter(!metaecosystem_type == "S_L") %>%
filter ( disturbance == "high") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = metaecosystem_type)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = high",
fill = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_colour_continuous(guide = "none") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large-large",
"medium-medium",
"small-small")) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional %>%
filter(disturbance == "low") %>%
filter(!metaecosystem_type == "S_L") %>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = interaction(day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(title = "Disturbance = low",
x = "Day",
y = "Local bioarea (something/μl)",
fill = "") +
#scale_fill_discrete(labels = c("isolated large", "isolated medium", "isolated small")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("large-large",
"medium-medium",
"small-small")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional %>%
filter(disturbance == "high") %>%
filter(!metaecosystem_type == "S_L") %>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = interaction(day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(title = "Disturbance = high",
x = "Day",
y = "Local bioarea (something/μl)",
fill = "") +
#scale_fill_discrete(labels = c("isolated large", "isolated medium", "isolated small")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("large-large",
"medium-medium",
"small-small")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
Interesting. It seems like there’s not much difference between the medium-medium and the large-large.
### --- SINGLE ECOSYSTEMS FOR ALL SMALL PATCHES --- ###
ds_biomass %>%
filter(disturbance == "low") %>%
#filter(eco_metaeco_type == "S (S_S)" | eco_metaeco_type == "S (S_L)") %>%
filter(patch_size == "S") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = eco_metaeco_type)) +
geom_line(stat = "summary", fun = "mean") +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = low",
linetype = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("small isolated",
"small connected to small",
"small connected to large")) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_biomass %>%
filter(disturbance == "high") %>%
#filter(eco_metaeco_type == "S (S_S)" | eco_metaeco_type == "S (S_L)") %>%
filter(patch_size == "S") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = eco_metaeco_type)) +
geom_line(stat = "summary", fun = "mean") +
labs(title = "Disturbance = high",
x = "Day",
y = "Local bioarea (something/μl)",
linetype = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("small isolated",
"small connected to small",
"small connected to large")) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
### --- BOXPLOTS FOR SMALL PATCHES CONNTED TO S OR L --- ###
local_small_low_plot = ds_biomass %>%
filter(disturbance == "low") %>%
filter(patch_size == "S") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day,eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = low",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("small isolated",
"small connected to small",
"small connected to large")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
local_small_low_plot
ds_biomass %>%
filter(disturbance == "high") %>%
filter(patch_size == "S") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day,eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(title = "Disturbance = high",
x = "Day",
y = "Local bioarea (something/μl)",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("small isolated",
"small connected to small",
"small connected to large")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
Let’s start from the largest model.
Should we
### --- SINGLE PATCHES --- ###
ds_biomass %>%
filter(disturbance == "low") %>%
filter(patch_size == "L") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = eco_metaeco_type)) +
geom_line(stat = "summary", fun = "mean") +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = low",
linetype = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large isolated",
"large connected to large",
"large connected to small"))
ds_biomass %>%
filter(disturbance == "high") %>%
filter(patch_size == "L") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = eco_metaeco_type)) +
geom_line(stat = "summary", fun = "mean") +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = high",
linetype = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large isolated",
"large connected to large",
"large connected to small"))
### --- BOXPLOTS --- ###
ds_biomass %>%
filter(disturbance == "low") %>%
filter(patch_size == "L") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day,eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = low",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("large isolated",
"large connected to large",
"large connected to small"))
local_large_low_plot = ds_biomass %>%
filter(disturbance == "high") %>%
filter(patch_size == "L") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day,eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = high",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("large isolated",
"large connected to large",
"large connected to small"))
Does M have an effect?
full_model = lmer(bioarea_per_volume ~
metaecosystem_type * disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L" |
eco_metaeco_type == "L (L_L)") %>%
filter(time_point >= 2),
REML = FALSE)
no_metaeco_type_model = lmer(bioarea_per_volume ~
disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L" |
eco_metaeco_type == "L (L_L)") %>%
filter(time_point >= 2),
REML = FALSE)
anova(full_model, no_metaeco_type_model)
## Data: ds_biomass %>% filter(eco_metaeco_type == "L" | eco_metaeco_type == ...
## Models:
## no_metaeco_type_model: bioarea_per_volume ~ disturbance + (1 | system_nr) + (1 | day)
## full_model: bioarea_per_volume ~ metaecosystem_type * disturbance + (1 | system_nr) + (1 | day)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_metaeco_type_model 5 3879.5 3896.9 -1934.7 3869.5
## full_model 7 3876.5 3900.9 -1931.2 3862.5 6.9701 2 0.03065
##
## no_metaeco_type_model
## full_model *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yes.
Does M have an effect?
full_model = lmer(bioarea_per_volume ~
metaecosystem_type * disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L" |
eco_metaeco_type == "L (S_L)") %>%
filter(time_point >= 2),
REML = FALSE)
no_metaeco_type_model = lmer(bioarea_per_volume ~
disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L" |
eco_metaeco_type == "L (S_L)") %>%
filter(time_point >= 2),
REML = FALSE)
anova(full_model, no_metaeco_type_model)
## Data: ds_biomass %>% filter(eco_metaeco_type == "L" | eco_metaeco_type == ...
## Models:
## no_metaeco_type_model: bioarea_per_volume ~ disturbance + (1 | system_nr) + (1 | day)
## full_model: bioarea_per_volume ~ metaecosystem_type * disturbance + (1 | system_nr) + (1 | day)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_metaeco_type_model 5 2566.9 2582.2 -1278.4 2556.9
## full_model 7 2559.7 2581.2 -1272.9 2545.7 11.151 2 0.00379
##
## no_metaeco_type_model
## full_model **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yes.
Does M have an effect?
full_model = lmer(bioarea_per_volume ~
metaecosystem_type * disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L (L_L)" |
eco_metaeco_type == "L (S_L)") %>%
filter(time_point >= 2),
REML = FALSE)
no_metaeco_type_model = lmer(bioarea_per_volume ~
disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L (L_L)" |
eco_metaeco_type == "L (S_L)") %>%
filter(time_point >= 2),
REML = FALSE)
anova(full_model, no_metaeco_type_model)
## Data: ds_biomass %>% filter(eco_metaeco_type == "L (L_L)" | eco_metaeco_type == ...
## Models:
## no_metaeco_type_model: bioarea_per_volume ~ disturbance + (1 | system_nr) + (1 | day)
## full_model: bioarea_per_volume ~ metaecosystem_type * disturbance + (1 | system_nr) + (1 | day)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_metaeco_type_model 5 3835.5 3852.9 -1912.8 3825.5
## full_model 7 3832.2 3856.6 -1909.1 3818.2 7.2508 2 0.02664
##
## no_metaeco_type_model
## full_model *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Yes.
ds_biomass %>%
filter ( disturbance == "low") %>%
filter(metaecosystem == "no") %>%
group_by (system_nr, day, patch_size) %>%
summarise(mean_bioarea_per_volume_across_videos = mean(bioarea_per_volume)) %>%
ggplot (aes(x = day,
y = mean_bioarea_per_volume_across_videos,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = patch_size)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = low",
fill = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_colour_continuous(guide = "none") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large isolated",
"medium isolated",
"small isolated"))
ds_biomass %>%
filter ( disturbance == "high") %>%
filter(metaecosystem == "no") %>%
group_by (system_nr, day, patch_size) %>%
summarise(mean_bioarea_per_volume_across_videos = mean(bioarea_per_volume)) %>%
ggplot (aes(x = day,
y = mean_bioarea_per_volume_across_videos,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = patch_size)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = low",
fill = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_colour_continuous(guide = "none") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large isolated",
"medium isolated",
"small isolated"))
ds_biomass %>%
filter(disturbance == "low") %>%
filter(metaecosystem == "no") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day, patch_size),
fill = patch_size)) +
geom_boxplot() +
labs(title = "Disturbance = low",
x = "Day",
y = "Local bioarea (something/μl)",
fill = "") +
scale_fill_discrete(labels = c("isolated large", "isolated medium", "isolated small")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6))
ds_biomass %>%
filter(disturbance == "high") %>%
filter(metaecosystem == "no") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day, patch_size),
fill = patch_size)) +
geom_boxplot() +
labs(title = "Disturbance = high",
x = "Day",
y = "Local bioarea (something/μl)",
fill = "") +
scale_fill_discrete(labels = c("isolated large", "isolated medium", "isolated small")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6))
We want to know if there was a systematic bias in the evaporation of different treatments (disturbance, patch size) and whether evaporation changed across time. My expectation would be that we would see a difference among the exchanges 2,3 and the exchanges 4,5,6. This is because in exchange 2,3 cultures were microwaved in 15 tubes for 3 minutes and in exchange 4,5,6 cultures were microwaved in 4 tubes for only 1 minute.
#Columns: exchange & evaporation
ds_for_evaporation = gather(ds_for_evaporation,
key = exchange,
value = evaporation,
water_add_after_t2:water_add_after_t6)
ds_for_evaporation[ds_for_evaporation == "water_add_after_t2"] = "2"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t3"] = "3"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t4"] = "4"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t5"] = "5"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t6"] = "6"
ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] = ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] / 2 #This is because exchange contained the topping up of two exchanges
ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] = ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] + 2 #We need to add 2 ml to the evaporation that happened at the exchange events 1 and 2. This is because we already added 1 ml of water at exchange 1 and 1 ml of water at exchange 2.
#Column: nr_of_tubes_in_rack
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 1] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 2] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 3] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 4] = 4
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 5] = 4
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 6] = 4
ds_for_evaporation %>%
filter(disturbance == disturbance) %>%
ggplot(aes(x = as.character(nr_of_tubes_in_rack),
y = evaporation)) +
geom_boxplot() +
labs(x = "Number of tubes in rack",
y = "Evaporation (ml)")
ds_for_evaporation %>%
filter(disturbance == disturbance) %>%
ggplot(aes(x = as.character(patch_size),
y = evaporation)) +
geom_boxplot() +
labs(x = "Patch size",
y = "Evaporation (ml)")
ds_for_evaporation %>%
filter(disturbance == disturbance) %>%
ggplot(aes(x = as.character(day),
y = evaporation)) +
geom_boxplot() +
labs(x = "Day",
y = "Evaporation (ml)")
ds_for_evaporation %>%
filter(disturbance == disturbance) %>%
ggplot(aes(x = disturbance,
y = evaporation)) +
geom_boxplot() +
labs(x = "Disturbance",
y = "Evaporation (ml)")
It seems like there is no real difference across time, disturbance, or patch type. However, we could also run a mixed effect model to show that they do not.
I’m not going to run it because it takes a long time.
mixed.model = lmer(evaporation ~
patch_size * disturbance * exchange +
(exchange | culture_ID),
data = ds_for_evaporation,
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
null.model = lm(evaporation ~
1,
data = ds_for_evaporation)
anova(mixed.model, null.model)
culture_info = read.csv(here("data", "PatchSizePilot_culture_info.csv"), header = TRUE)
load(here("data", "morphology", "t0.RData"));t0 = morph_mvt
load(here("data", "morphology", "t1.RData"));t1 = morph_mvt
load(here("data", "morphology", "t2.RData"));t2 = morph_mvt
load(here("data", "morphology", "t3.RData"));t3 = morph_mvt
load(here("data", "morphology", "t4.RData"));t4 = morph_mvt
load(here("data", "morphology", "t5.RData"));t5 = morph_mvt
load(here("data", "morphology", "t6.RData"));t6 = morph_mvt
load(here("data", "morphology", "t7.RData"));t7 = morph_mvt
rm(morph_mvt)
### --- Tidy t0 - t7 data-sets --- ###
#Column: time
t0$time = NA
t1$time = NA
#Column: replicate_video
t0$replicate_video[t0$file == "sample_00001"] = 1
t0$replicate_video[t0$file == "sample_00002"] = 2
t0$replicate_video[t0$file == "sample_00003"] = 3
t0$replicate_video[t0$file == "sample_00004"] = 4
t0$replicate_video[t0$file == "sample_00005"] = 5
t0$replicate_video[t0$file == "sample_00006"] = 6
t0$replicate_video[t0$file == "sample_00007"] = 7
t0$replicate_video[t0$file == "sample_00008"] = 8
t0$replicate_video[t0$file == "sample_00009"] = 9
t0$replicate_video[t0$file == "sample_00010"] = 10
t0$replicate_video[t0$file == "sample_00011"] = 11
t0$replicate_video[t0$file == "sample_00012"] = 12
t1$replicate_video = 1 #In t1 I took only 1 video/culture
t2$replicate_video = 1 #In t2 I took only 1 video/culture
t3$replicate_video = 1 #In t3 I took only 1 video/culture
t4$replicate_video = 1 #In t4 I took only 1 video/culture
t5$replicate_video = 1 #In t5 I took only 1 video/culture
t6 = t6 %>% rename(replicate_video = replicate)
t7 = t7 %>% rename(replicate_video = replicate)
### --- Create ds_body_size dataset --- ###
long_t0 = t0 %>% slice(rep(1:n(), max(culture_info$culture_ID)))
ID_vector = NULL
ID_vector_elongating = NULL
for (ID in 1:max(culture_info$culture_ID)){
ID_vector = rep(ID, times = nrow(t0))
ID_vector_elongating = c(ID_vector_elongating, ID_vector)
}
long_t0$culture_ID = ID_vector_elongating
t0 = merge(culture_info,long_t0, by="culture_ID"); rm(long_t0)
t1 = merge(culture_info,t1,by="culture_ID")
t2 = merge(culture_info,t2,by="culture_ID")
t3 = merge(culture_info,t3,by="culture_ID")
t4 = merge(culture_info,t4,by="culture_ID")
t5 = merge(culture_info,t5,by="culture_ID")
t6 = merge(culture_info,t6,by="culture_ID")
t7 = merge(culture_info,t7,by="culture_ID")
ds_body_size = rbind(t0, t1, t2, t3, t4, t5, t6, t7); rm(t0, t1, t2, t3, t4, t5, t6, t7)
### --- Tidy ds_body_size data-set --- ###
#Column: day
ds_body_size$day = ds_body_size$time_point;
ds_body_size$day[ds_body_size$day=="t0"] = "0"
ds_body_size$day[ds_body_size$day=="t1"] = "4"
ds_body_size$day[ds_body_size$day=="t2"] = "8"
ds_body_size$day[ds_body_size$day=="t3"] = "12"
ds_body_size$day[ds_body_size$day=="t4"] = "16"
ds_body_size$day[ds_body_size$day=="t5"] = "20"
ds_body_size$day[ds_body_size$day=="t6"] = "24"
ds_body_size$day[ds_body_size$day=="t7"] = "28"
ds_body_size$day = as.numeric(ds_body_size$day)
#Column: time point
ds_body_size$time_point[ds_body_size$time_point=="t0"] = 0
ds_body_size$time_point[ds_body_size$time_point=="t1"] = 1
ds_body_size$time_point[ds_body_size$time_point=="t2"] = 2
ds_body_size$time_point[ds_body_size$time_point=="t3"] = 3
ds_body_size$time_point[ds_body_size$time_point=="t4"] = 4
ds_body_size$time_point[ds_body_size$time_point=="t5"] = 5
ds_body_size$time_point[ds_body_size$time_point=="t6"] = 6
ds_body_size$time_point[ds_body_size$time_point=="t7"] = 7
ds_body_size$time_point = as.character(ds_body_size$time_point)
#Column: eco_metaeco_type
ds_body_size$eco_metaeco_type = factor(ds_body_size$eco_metaeco_type,
levels=c('S', 'S (S_S)', 'S (S_L)', 'M', 'M (M_M)', 'L', 'L (L_L)', 'L (S_L)'))
#Select useful columns
ds_body_size = ds_body_size %>%
select(culture_ID,
patch_size,
disturbance,
metaecosystem_type,
mean_area,
replicate_video,
day,
metaecosystem,
system_nr,
eco_metaeco_type)
#Reorder columns
ds_body_size = ds_body_size[, c("culture_ID",
"system_nr",
"disturbance",
"day",
"patch_size",
"metaecosystem",
"metaecosystem_type",
"eco_metaeco_type",
"replicate_video",
"mean_area")]
As in Jacquet, Gounand, and Altermatt (2020) I will create 12 size classes.
#### --- PARAMETERS & INITIALISATION --- ###
nr_of_size_classes = 12
largest_size = max(ds_body_size$mean_area)
size_class_width = largest_size/nr_of_size_classes
size_class = NULL
### --- CREATE DATASET --- ###
size_class_boundaries = seq(0, largest_size, by = size_class_width)
for (class in 1:nr_of_size_classes){
bin_lower_limit = size_class_boundaries[class]
bin_upper_limit = size_class_boundaries[class+1]
size_input = (size_class_boundaries[class] + size_class_boundaries[class + 1])/2
size_class[[class]] = ds_body_size%>%
filter(bin_lower_limit <= mean_area) %>%
filter(mean_area <= bin_upper_limit) %>%
group_by(culture_ID,
system_nr,
disturbance,
day,
patch_size,
metaecosystem,
metaecosystem_type,
eco_metaeco_type,
replicate_video) %>% #Group by video
summarise(mean_abundance_across_videos = n()) %>%
group_by(culture_ID,
system_nr,
disturbance,
day,
patch_size,
metaecosystem,
metaecosystem_type,
eco_metaeco_type) %>% #Group by ID
summarise(abundance = mean(mean_abundance_across_videos)) %>%
mutate(log_abundance = log(abundance)) %>%
mutate(size_class = class) %>%
mutate(size = size_input) %>%
mutate(log_size = log(size))
}
ds_classes = rbind(size_class[[1]], size_class[[2]], size_class[[3]], size_class[[4]],
size_class[[5]], size_class[[6]], size_class[[7]], size_class[[8]],
size_class[[9]], size_class[[10]], size_class[[11]], size_class[[12]],)
The code for plotting the size class distribution across small patches doesn’t work for some reason. It displays all plots as if they were the same.

Regional bioarea density (mean bioarea density between two patches) in meta-ecosystems of the same total area, but whose two patches have either the same size (red, medium-medium meta-ecosystem) or that have a smaller and larger patch (blue, small-large meta-ecosystem). Points represent the mean, error bars represent the standard deviation. For clarity, only the low disturbance treatment is shown here. See the Appendix for equivalent the figure of the high disturbance treatment.
Local biomass production in patches that are either connnected to a patch of the same size (green, small patches of small-small meta-ecosystemes) or to a patch of larger size (orange, small patches of small-large meta-ecosystems). Points represent the mean, error bars represent the standard deviation. For clarity, only the low disturbance treatment is shown here. See the Appendix for equivalent the figure of the high disturbance treatment. Problem: the small patches in green actually have different sample size than the one in orange. Does it still make sense to include the standard deviation?
evaporation.test = read.csv(here("data", "evaporation_test","evaporation_test_right.csv"), header = TRUE)
evaporation.test %>%
ggplot(aes (x = as.character(water_pipetted),
y = weight_water_evaporated,
group = interaction(water_pipetted, as.character(rack)),
fill = as.character(rack))) +
geom_boxplot() +
labs(x = "Water volume (ml)" ,
y = "Evaporation (g)",
fill = "Rack replicate")
evaporation.test = read.csv(here("data", "evaporation_test", "evaporation_test_fill_nofill.csv"), header = TRUE)
evaporation.test %>%
ggplot(aes (x = all_tubes_water,
y = weight_water_evaporated)) +
geom_boxplot() +
labs(x = "Water in the other 10 tubes" ,
y = "Evaporation (g)",
caption = "When all tubes were filled, they were filled with 6.75 ml of deionised water (I think, but I need to check in my lab book.")
## Time difference of 1.111363 mins
To build the mixed effect models we will use the R package lme4. See page 6 of this PDF to know more about the syntaxis of this package.
To do model diagnostics of mixed effect models, I’m going to look at the following two plots (as suggested by Zuur et al. (2009), page 487):
Quantile-quantile plots
Partial residual plots
The effect size of the explaining variables is calculated in the
mixed effect models as marginal and conditional r squared. The marginal
r squared is how much variance is explained by the fixed effects. The
conditional r squared is how much variance is explained by the fixed and
the random effects. The marginal and conditional r squared are
calculated using the package MuMIn. The computation is
based on the methods of Nakagawa, Johnson, and
Schielzeth (2017). For the coding and interpretation of these r
squared check the documentation
for the r.squaredGLMM function
See for the interaction syntaxis this link.
I am starting from the largest model and then simplifying because … (see statistical modelling course at ETH).
I am going to select the best model according to AIC. BIC is better for understanding and AIC for predicting. Halsey (2019) also suggests this approach instead of p values. I’m going to use AIC because I’m interested in knowing how much meta-ecosystem type contributed to the overall regional biomass.